home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / derivsig.pro < prev    next >
Text File  |  1997-07-08  |  4KB  |  106 lines
  1. ; $Id: derivsig.pro,v 1.3 1996/12/17 23:43:58 lubos Exp $
  2. Function Derivsig, X, Y, sigx, sigy
  3. ;+
  4. ; NAME:
  5. ;    DERIVSIG
  6. ;
  7. ; PURPOSE:
  8. ;    This function computes the standard deviation of a derivative
  9. ;       as found by the DERIV function, using the input variables of
  10. ;    DERIV and the standard deviations of those input variables.
  11. ;
  12. ; CATEGORY:
  13. ;    Numerical analysis.
  14. ;
  15. ; CALLING SEQUENCE:
  16. ;    sigDy = Derivsig(sigy)        ;sigma(Dy(i)/di), point spacing = 1.
  17. ;    sigDy = Derivsig(X,Y,sigx,sigy) ;sigma(Dy/Dx), unequal point spacing.
  18. ;
  19. ; INPUTS:
  20. ;    Y:    The variable to be differentiated. Omit if X is omitted.
  21. ;    X:    The Variable to differentiate with respect to. If omitted,
  22. ;        unit spacing is assumed for Y, i.e. X(i) = i.
  23. ;       sigy:    The standard deviation of Y. (Vector if used alone in
  24. ;        call, vector or constant if used with other parameters)
  25. ;       sigx:    The standard deviation of X (either vector or constant).
  26. ;        Use "0.0" if the abscissa is exact; omit if X is omitted.
  27. ;
  28. ; OPTIONAL INPUT PARAMETERS:
  29. ;    As above.
  30. ;
  31. ; OUTPUTS:
  32. ;    This function returns the standard deviation of the derivative.
  33. ;
  34. ; COMMON BLOCKS:
  35. ;    None.
  36. ;
  37. ; SIDE EFFECTS:
  38. ;    None.
  39. ;
  40. ; RESTRICTIONS:
  41. ;    None.
  42. ;
  43. ; PROCEDURE:
  44. ;    See Bevington, "Data Analysis and Reduction for the Physical
  45. ;           Sciences," McGraw-Hill (1969), Chap 4.
  46. ;
  47. ; MODIFICATION HISTORY:
  48. ;       Written by Richard Bonomo at the University of Wisconsin - Madison
  49. ;       department of Electrical and Computer Engineering, July, 1991.
  50. ;    "DERIV" written by DMS, Aug, 1984.
  51. ;-
  52. ;
  53.  
  54. on_error,2              ;Return to caller if an error occurs
  55. prms=n_params(0)
  56. n = n_elements(x)
  57. if (n lt 3) and (prms gt 1) then message, 'X must have at least 3 points'
  58. if (n lt 3) and (prms eq 1) then message, $
  59.     'sigy must be a vector of at least 3 points if used alone'
  60.  
  61. if ((prms ne 1) and (prms ne 4)) then message,$
  62.    'function DERIVSIG must be called with either 1 or 4 parameters'
  63.  
  64. if prms eq 1 then begin ; unit spacing assumed
  65.         sigy=x
  66.         if n_elements(sigy) eq 1 then sigy=fltarr(n) + sigy
  67.         sigd=sqrt(0.25*(shift(sigy,-1)*shift(sigy,-1) + $
  68.          shift(sigy,1)*shift(sigy,1)))
  69.         sigd[0]=sqrt(0.25*(sigy[0]^2*9.0 + sigy[1]^2*16.0 + sigy[2]^2))
  70.         sigd[n-1]=sqrt(0.25*(sigy[n-1]^2*9.0 + sigy[n-2]*16.0 + sigy[n-3]))
  71. endif
  72. if prms eq 4 then begin
  73.         if n ne n_elements(y) then message,'Vectors must have same size'
  74.         if n_elements(sigy) eq 1 then sigy=fltarr(n) + sigy
  75.         nix=n_elements(sigx)
  76.         if (nix eq 1) and (sigx[0] ne 0.0) then sigx=fltarr(n) + sigx
  77.         nix=n_elements(sigx)
  78.         if (nix ne n) and (nix ne 1) then message,$
  79.            'sigx vector must have the same length as X, or be a scalar'
  80.         dsq=shift(x,-1)-shift(x,1)
  81.         dsq=dsq*dsq
  82.         dy=shift(y,-1)-shift(y,1)
  83.         sigd=(shift(sigy,-1)*shift(sigy,-1) + shift(sigy,1)*shift(sigy,1))/dsq
  84.         if (nix ne 1) then sigd=sigd + (shift(sigx,-1)^2*dy^2 + $
  85.          shift(sigx,1)^2*dy^2)/(dsq*dsq)
  86.         sigd=sqrt(sigd)
  87.         dsq=x[2]-x[0]
  88.         dsq=dsq*dsq
  89. ;
  90.         sigd[0]=(sigy[0]^2*9.0 + sigy[1]^2*16.0 + sigy[2]^2)/dsq
  91.         if (nix ne 1) then sigd[0] = sigd[0] + (sigx[2]^2*(3.0*y[0] - $
  92.          4.0*y[1] + y[2])^2 + sigx[0]^2*(-3.0*y[0] + 4.0*y[1] - y[2])^2)/$
  93.          (dsq*dsq)
  94.         sigd[0]=sqrt(sigd[0])
  95. ;
  96.         dsq=x[n-1]-x[n-3]
  97.         dsq=dsq*dsq
  98.         sigd[n-1]=(sigy[n-1]^2*9.0 + sigy[n-2]^2*16.0 + sigy[n-3]^2)/dsq
  99.         if (nix ne 1) then sigd[n-1] = sigd[n-1] + (sigx[n-1]^2*(-3.0*y[n-1]$
  100.          + 4.0*y[n-2] - y[n-3])^2 + sigx[n-3]^2*(3.0*y[n-1] - 4.0*y[n-2]$
  101.          + y[n-3])^2)/(dsq*dsq)
  102.         sigd[n-1]=sqrt(sigd[n-1])
  103. endif
  104. return, sigd
  105. end
  106.